
# import packages and dataset
import seaborn as sns
%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
url = 'https://raw.githubusercontent.com/christianclayton7/ibm-attrition-prediction/main/ibm_attrition_data.csv'
df = pd.read_csv(url)
pd.options.display.max_columns = None


# preview data
df.head()


# look at columns
df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1470 entries, 0 to 1469
Data columns (total 35 columns):
 #   Column                    Non-Null Count  Dtype 
---  ------                    --------------  ----- 
 0   Age                       1470 non-null   int64 
 1   Attrition                 1470 non-null   object
 2   BusinessTravel            1470 non-null   object
 3   DailyRate                 1470 non-null   int64 
 4   Department                1470 non-null   object
 5   DistanceFromHome          1470 non-null   int64 
 6   Education                 1470 non-null   int64 
 7   EducationField            1470 non-null   object
 8   EmployeeCount             1470 non-null   int64 
 9   EmployeeNumber            1470 non-null   int64 
 10  EnvironmentSatisfaction   1470 non-null   int64 
 11  Gender                    1470 non-null   object
 12  HourlyRate                1470 non-null   int64 
 13  JobInvolvement            1470 non-null   int64 
 14  JobLevel                  1470 non-null   int64 
 15  JobRole                   1470 non-null   object
 16  JobSatisfaction           1470 non-null   int64 
 17  MaritalStatus             1470 non-null   object
 18  MonthlyIncome             1470 non-null   int64 
 19  MonthlyRate               1470 non-null   int64 
 20  NumCompaniesWorked        1470 non-null   int64 
 21  Over18                    1470 non-null   object
 22  OverTime                  1470 non-null   object
 23  PercentSalaryHike         1470 non-null   int64 
 24  PerformanceRating         1470 non-null   int64 
 25  RelationshipSatisfaction  1470 non-null   int64 
 26  StandardHours             1470 non-null   int64 
 27  StockOptionLevel          1470 non-null   int64 
 28  TotalWorkingYears         1470 non-null   int64 
 29  TrainingTimesLastYear     1470 non-null   int64 
 30  WorkLifeBalance           1470 non-null   int64 
 31  YearsAtCompany            1470 non-null   int64 
 32  YearsInCurrentRole        1470 non-null   int64 
 33  YearsSinceLastPromotion   1470 non-null   int64 
 34  YearsWithCurrManager      1470 non-null   int64 
dtypes: int64(26), object(9)
memory usage: 402.1+ KB


# look at summary stats of numeric features
pd.options.display.precision = 2
df.describe()


# look at proportion of each level in each category feature
for col in df.select_dtypes(include='object').columns: 
    print("\n", df[col].value_counts(normalize = True).round(2))

 Attrition
No     0.84
Yes    0.16
Name: proportion, dtype: float64

 BusinessTravel
Travel_Rarely        0.71
Travel_Frequently    0.19
Non-Travel           0.10
Name: proportion, dtype: float64

 Department
Research & Development    0.65
Sales                     0.30
Human Resources           0.04
Name: proportion, dtype: float64

 EducationField
Life Sciences       0.41
Medical             0.32
Marketing           0.11
Technical Degree    0.09
Other               0.06
Human Resources     0.02
Name: proportion, dtype: float64

 Gender
Male      0.6
Female    0.4
Name: proportion, dtype: float64

 JobRole
Sales Executive              0.22
Research Scientist           0.20
Laboratory Technician        0.18
Manufacturing Director       0.10
Healthcare Representative    0.09
Manager                      0.07
Sales Representative         0.06
Research Director            0.05
Human Resources              0.04
Name: proportion, dtype: float64

 MaritalStatus
Married     0.46
Single      0.32
Divorced    0.22
Name: proportion, dtype: float64

 Over18
Y    1.0
Name: proportion, dtype: float64

 OverTime
No     0.72
Yes    0.28
Name: proportion, dtype: float64


# view distribution of numeric features (excepting ones with very few unique values) by attrition
drop_cols = ['StandardHours','EmployeeCount','EmployeeNumber',
             'WorkLifeBalance','RelationshipSatisfaction','JobSatisfaction','JobInvolvement','EnvironmentSatisfaction',
            'PerformanceRating', 'StockOptionLevel']
num_cols = df.drop(drop_cols, axis = 1).select_dtypes(include='int64').columns
for col in num_cols:
    sns.catplot(kind = 'boxen', data = df, x = col, y = 'Attrition')
    plt.show()

/Users/christianclayton/anaconda3/lib/python3.11/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
  self._figure.tight_layout(*args, **kwargs)

/Users/christianclayton/anaconda3/lib/python3.11/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
  self._figure.tight_layout(*args, **kwargs)

/Users/christianclayton/anaconda3/lib/python3.11/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
  self._figure.tight_layout(*args, **kwargs)

/Users/christianclayton/anaconda3/lib/python3.11/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
  self._figure.tight_layout(*args, **kwargs)

/Users/christianclayton/anaconda3/lib/python3.11/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
  self._figure.tight_layout(*args, **kwargs)

/Users/christianclayton/anaconda3/lib/python3.11/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
  self._figure.tight_layout(*args, **kwargs)

/Users/christianclayton/anaconda3/lib/python3.11/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
  self._figure.tight_layout(*args, **kwargs)

/Users/christianclayton/anaconda3/lib/python3.11/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
  self._figure.tight_layout(*args, **kwargs)

/Users/christianclayton/anaconda3/lib/python3.11/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
  self._figure.tight_layout(*args, **kwargs)

/Users/christianclayton/anaconda3/lib/python3.11/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
  self._figure.tight_layout(*args, **kwargs)


# view proportion of leaves in each level of ordinal variable (numeric with few unique values)
cat_cols = df[['WorkLifeBalance','RelationshipSatisfaction','JobSatisfaction','JobInvolvement','EnvironmentSatisfaction',
            'PerformanceRating', 'StockOptionLevel']].columns
for col in cat_cols:
    temp = df.groupby(col).value_counts(['Attrition'], normalize = True).reset_index()
    sns.barplot(data = temp,
            x = col,
            y = 'proportion',
            hue = 'Attrition')
    plt.show()


# view proportion of leaves in each category level of category variables
cat_cols = df.drop('Over18', axis = 1).select_dtypes(include='object').columns
cat_cols = cat_cols[cat_cols != 'Attrition']
for col in cat_cols:
    temp = df.groupby(col).value_counts(['Attrition'], normalize = True).reset_index()
    sns.barplot(data = temp,
            x = col,
            y = 'proportion',
            hue = 'Attrition')
    plt.xticks(rotation=90)
    plt.show()


# function for sorting the correlation matrix
import scipy
import scipy.cluster.hierarchy as sch
def cluster_corr(corr_array, inplace=False):
    """
    Rearranges the correlation matrix, corr_array, so that groups of highly 
    correlated variables are next to eachother 
    
    Parameters
    ----------
    corr_array : pandas.DataFrame or numpy.ndarray
        a NxN correlation matrix 
        
    Returns
    -------
    pandas.DataFrame or numpy.ndarray
        a NxN correlation matrix with the columns and rows rearranged
    """
    pairwise_distances = sch.distance.pdist(corr_array)
    linkage = sch.linkage(pairwise_distances, method='complete')
    cluster_distance_threshold = pairwise_distances.max()/2
    idx_to_cluster_array = sch.fcluster(linkage, cluster_distance_threshold, 
                                        criterion='distance')
    idx = np.argsort(idx_to_cluster_array)
    
    if not inplace:
        corr_array = corr_array.copy()
    
    if isinstance(corr_array, pd.DataFrame):
        return corr_array.iloc[idx, :].T.iloc[idx, :]
    return corr_array[idx, :][:, idx]

# get correlations and sort them
corr = cluster_corr(df.drop(['Over18','StandardHours','EmployeeCount','EmployeeNumber'],axis = 1).select_dtypes(include='int64').corr())

# plot correlations in heatmap
sns.heatmap(corr)

<Axes: >


# check for number of significant principle components on year count/experience features
from sklearn.decomposition import PCA
pca = PCA()
experience_pca = pca.fit_transform(df[['Age','JobLevel', 'TotalWorkingYears','YearsAtCompany', 'YearsInCurrentRole',
                                'YearsSinceLastPromotion','YearsWithCurrManager']])

def cum_pca_var(n_comps, cum_exp_var_ratio):
    plt.plot(range(1,n_comps + 1), cum_exp_var_ratio)
    plt.annotate('cumulative', (n_comps + 1,1), (n_comps+ 1,.9))
    plt.xticks(range(n_comps + 1))
    plt.ylabel('% of variance explained')
    plt.xlabel('Principle component')
    plt.title('Cumulative % of variance explained by additional components')
    plt.show()

cum_pca_var(pca.n_components_, 
            pca.explained_variance_ratio_.cumsum()
           )


# check for number of significant principle components on compensation features
pca = PCA()
sent_pca = pca.fit_transform(df[['DailyRate','HourlyRate','MonthlyIncome', 'MonthlyRate']])
cum_pca_var(pca.n_components_, 
            pca.explained_variance_ratio_.cumsum()
           )


# check for number of significant principle components on sentiment features
pca = PCA()
sent_pca = pca.fit_transform(df[['WorkLifeBalance','RelationshipSatisfaction','JobSatisfaction',
                                 'JobInvolvement','EnvironmentSatisfaction']])
cum_pca_var(pca.n_components_, 
            pca.explained_variance_ratio_.cumsum()
           )


# get 2 principle components of the experience features
pca = PCA(n_components = 2)
experience_pca = pca.fit_transform(df[['Age','JobLevel', 'TotalWorkingYears','YearsAtCompany', 'YearsInCurrentRole',
                                'YearsSinceLastPromotion','YearsWithCurrManager']])
experience_pca = pd.DataFrame(experience_pca)
experience_pca.columns = ['Experience1', 'Experience2']

# get 2 principle components of the compensation features
pca = PCA(n_components = 2)
compensation_pca = pca.fit_transform(df[['DailyRate','HourlyRate','MonthlyIncome', 'MonthlyRate']])
compensation_pca = pd.DataFrame(compensation_pca)
compensation_pca.columns = ['Compensation1', 'Compensation2']

# get 3 principle components of the sentiment features
pca = PCA(n_components = 3)
sentiment_pca = pca.fit_transform(df[['WorkLifeBalance','RelationshipSatisfaction','JobSatisfaction',
                                 'JobInvolvement','EnvironmentSatisfaction']])
sentiment_pca = pd.DataFrame(sentiment_pca)
sentiment_pca.columns = ['Sentiment1', 'Sentiment2', 'Sentiment3']


# turn categorical features into 1's and 0's for sklearn to be able to interpret and use
df_cats = pd.get_dummies(df.drop(['Over18'], axis = 1).select_dtypes(include='object'), drop_first = True)


# combine pca features, categorical dummy features, and remaining useful numeric features
drop_cols = ['StandardHours','EmployeeCount','EmployeeNumber', 'Over18', 
             'WorkLifeBalance','RelationshipSatisfaction','JobSatisfaction','JobInvolvement','EnvironmentSatisfaction',
             'DailyRate','HourlyRate','MonthlyIncome', 'MonthlyRate',
             'Age','JobLevel', 'TotalWorkingYears','YearsAtCompany', 'YearsInCurrentRole',
             'YearsSinceLastPromotion','YearsWithCurrManager',
             'Attrition',
             'BusinessTravel', 'Department', 'EducationField', 'Gender', 'JobRole', 'MaritalStatus', 'OverTime']

# combine dataframes
df2 = pd.concat([df.drop(drop_cols, axis = 1),
          df_cats,
          experience_pca,
          compensation_pca,
          sentiment_pca],
         axis = 1)


# split into training and test sets
from sklearn.model_selection import train_test_split
X = df2.drop('Attrition_Yes', axis = 1).values
y = df2['Attrition_Yes'].values
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state = 7, test_size = .3, stratify = y)

# standardize features
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.fit_transform(X_test)


# get model packages
from sklearn.ensemble import GradientBoostingClassifier, RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import cross_val_score, KFold

# set up dictionary of model types to estimate
models = {'KNN':KNeighborsClassifier(n_neighbors = 5),
          'Logit':LogisticRegression(),
          'RF':RandomForestClassifier(n_estimators = 80, min_samples_leaf = 5, random_state = 7),
          'GBM':GradientBoostingClassifier(n_estimators = 80, max_depth = 3, subsample = .8, max_features = .2, random_state = 7)}

# set up empty list to catch CV training scores
results_accuracy = []
results_recall = []

# loop through CV scoring on each model type
for model in models.values():
    kf = KFold(n_splits = 6, random_state = 7, shuffle = True)
    cv_results_accuracy = cross_val_score(model, X_train_scaled, y_train, cv = kf)
    cv_results_recall = cross_val_score(model, X_train_scaled, y_train, cv = kf, scoring = 'recall')
    results_accuracy.append(cv_results_accuracy)
    results_recall.append(cv_results_recall)

# plot cv scores
plt.boxplot(results_accuracy, labels = models.keys())
plt.title('Total CV Accuracy')
plt.show()

# plot cv scores
plt.boxplot(results_recall, labels = models.keys())
plt.title('Total CV Recall (Accuracy on actual leaves)')
plt.show()


from sklearn.model_selection import GridSearchCV
kf = KFold(n_splits = 5, shuffle = True, random_state = 7)
param_grid = {'penalty' : ['l1', 'l2'],
              'C' : np.logspace(-4, 4, 5),
              'solver' : ['liblinear']}
logit = LogisticRegression()
logit_cv = GridSearchCV(logit, param_grid, cv = kf, scoring = 'recall')
logit_cv.fit(X_train_scaled, y_train)
print('top tuned parameters:', logit_cv.best_params_)
print('top tuned recall score:', round(logit_cv.best_score_,2))

top tuned parameters: {'C': 0.0001, 'penalty': 'l2', 'solver': 'liblinear'}
top tuned recall score: 0.4


from sklearn.metrics import recall_score
logit_best = LogisticRegression(penalty = 'l2', C = .0001, solver = 'liblinear')
logit_best.fit(X_train_scaled, y_train)
y_pred = logit_best.predict(X_test_scaled)
print('test data recall score:', round(recall_score(y_test, y_pred), 2))
print('test data overall accuracy', round(logit_best.score(X_test_scaled, y_test),2))

test data recall score: 0.38
test data overall accuracy 0.86


# inspect regression coefficients
coefs_df = pd.DataFrame({'Feature':list(df2.drop('Attrition_Yes', axis = 1).columns),
            'Coef':logit_best.coef_[0]}).sort_values('Coef',ascending = False)
coefs_df['Coef (abs)'] = coefs_df['Coef'].abs()
coefs_df = coefs_df.sort_values('Coef (abs)', ascending = False)
coefs_df['Attrition Odds Change'] = np.where(coefs_df['Coef'] < 0, 'Decrease','Increase')
condition = coefs_df['Feature'].isin(['OverTime_Yes', 'Compensation2', 'Compensation1', 'StockOptionLevel', 'Sentiment2',
                'Sentiment3', 'Sentiment1', 'BusinessTravel_Travel_Frequently', 'DistanceFromHome',
                'TrainingTimesLastYear', 'BusinessTravel_Travel_Rarely', 'PercentSalaryHike'])
coefs_df.loc[np.invert(condition), 'influence'] = 'cannot influence'
coefs_df.loc[condition, 'influence'] = 'can influence'



fig, (ax1, ax2) = plt.subplots(1,2)
plt.rcParams['figure.figsize'] = [12, 8]
sns.barplot(data = coefs_df[coefs_df['influence'] == 'can influence'], 
            x = 'Coef (abs)', 
            y = 'Feature', 
            hue = 'Attrition Odds Change', 
            dodge = False, 
            palette = 'RdBu',
           ax = ax1)
sns.barplot(data = coefs_df[coefs_df['influence'] == 'cannot influence'], 
            x = 'Coef (abs)', 
            y = 'Feature', 
            hue = 'Attrition Odds Change', 
            dodge = False, 
            palette = 'RdBu',
           ax = ax2)
fig.suptitle('Feature relationships to attrition odds by whether they can be influenced')
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
ax1.title.set_text('can influence')
ax2.title.set_text('cannot influence')
plt.show()


import plotly.io as pio
pio.renderers.default = 'notebook'


# get predicted attrition probabilities
y_pred_proba = pd.concat(
    [pd.DataFrame(logit_best.predict_proba(X_test_scaled)), 
     pd.DataFrame(logit_best.predict_proba(X_train_scaled))]
)

# get original columns
X = df.drop('Attrition', axis = 1).values
y = df['Attrition'].values
X_train_orig, X_test_orig, y_train_orig, y_test_orig = train_test_split(X, y, random_state = 7, test_size = .3, stratify = y)

# set up new dataframe with original columns and add predicted quit probabilities
pred_df = pd.concat(
    [pd.DataFrame(X_test_orig),
    pd.DataFrame(X_train_orig)]
)
pred_df = pd.concat([pd.DataFrame(y_pred_proba), pred_df], axis = 1)
pred_df.columns = ['Stay Probability', 'Leave Probability'] + list(df.drop('Attrition', axis = 1).columns)
pred_df['Prediction'] = np.where(pred_df['Leave Probability'] >= .5, 'Leave', 'Stay')
pred_df['YearlySalary'] = pred_df['MonthlyIncome'] * 12
pred_df['Leave Probability'] = pred_df['Leave Probability']

# visualize probleave against hard to replace characteristics
import plotly.express as px
#import plotly.io as pio
#pio.renderers.default = 'notebook'
fig = px.scatter(pred_df, 
                 x='Leave Probability', 
                 y="JobRole", 
                 color='Leave Probability', 
                 color_continuous_scale='Bluered',
                 facet_col="PerformanceRating",
                 facet_row='JobLevel',
                 width=1000, height=1400,
                 category_orders={"JobLevel": [5,4,3,2,1]},
                 title = 'Identifying hard-to-replace employees with high leave probability',
                 hover_data = {'PerformanceRating':False, 'JobRole':False, 'JobLevel':False,
                              'EmployeeNumber':True, 'Leave Probability':':.3f', 'Prediction':True, 'OverTime':True, 
                               'YearlySalary':True,
                             'StockOptionLevel':True, 'BusinessTravel':True,'DistanceFromHome':True})
fig.add_vline(x = .5, annotation_text="prediction boundary", annotation_position="bottom right" ,
              line_dash = 'dot')
fig.update_layout(coloraxis_showscale=False)
fig.show()

	Age	DailyRate	DistanceFromHome	Education	EmployeeCount	EmployeeNumber	EnvironmentSatisfaction	HourlyRate	JobInvolvement	JobLevel	JobSatisfaction	MonthlyIncome	MonthlyRate	NumCompaniesWorked	PercentSalaryHike	PerformanceRating	RelationshipSatisfaction	StandardHours	StockOptionLevel	TotalWorkingYears	TrainingTimesLastYear	WorkLifeBalance	YearsAtCompany	YearsInCurrentRole	YearsSinceLastPromotion	YearsWithCurrManager
count	1470.00	1470.00	1470.00	1470.00	1470.0	1470.00	1470.00	1470.00	1470.00	1470.00	1470.00	1470.00	1470.00	1470.00	1470.00	1470.00	1470.00	1470.0	1470.00	1470.00	1470.00	1470.00	1470.00	1470.00	1470.00	1470.00
mean	36.92	802.49	9.19	2.91	1.0	1024.87	2.72	65.89	2.73	2.06	2.73	6502.93	14313.10	2.69	15.21	3.15	2.71	80.0	0.79	11.28	2.80	2.76	7.01	4.23	2.19	4.12
std	9.14	403.51	8.11	1.02	0.0	602.02	1.09	20.33	0.71	1.11	1.10	4707.96	7117.79	2.50	3.66	0.36	1.08	0.0	0.85	7.78	1.29	0.71	6.13	3.62	3.22	3.57
min	18.00	102.00	1.00	1.00	1.0	1.00	1.00	30.00	1.00	1.00	1.00	1009.00	2094.00	0.00	11.00	3.00	1.00	80.0	0.00	0.00	0.00	1.00	0.00	0.00	0.00	0.00
25%	30.00	465.00	2.00	2.00	1.0	491.25	2.00	48.00	2.00	1.00	2.00	2911.00	8047.00	1.00	12.00	3.00	2.00	80.0	0.00	6.00	2.00	2.00	3.00	2.00	0.00	2.00
50%	36.00	802.00	7.00	3.00	1.0	1020.50	3.00	66.00	3.00	2.00	3.00	4919.00	14235.50	2.00	14.00	3.00	3.00	80.0	1.00	10.00	3.00	3.00	5.00	3.00	1.00	3.00
75%	43.00	1157.00	14.00	4.00	1.0	1555.75	4.00	83.75	3.00	3.00	4.00	8379.00	20461.50	4.00	18.00	3.00	4.00	80.0	1.00	15.00	3.00	3.00	9.00	7.00	3.00	7.00
max	60.00	1499.00	29.00	5.00	1.0	2068.00	4.00	100.00	4.00	5.00	4.00	19999.00	26999.00	9.00	25.00	4.00	4.00	80.0	3.00	40.00	6.00	4.00	40.00	18.00	15.00	17.00

Predicting employee attrition¶

Dataset¶

Validation:¶

Exploratory analysis¶

Principle Component Analysis¶

Reducing redundancy and noise between model features by consolidating information from many correlated features to a few uncorrelated features:¶

Feature Engineering, Transformations, and Train/Test Split¶

Model Fitting and Evaluation¶

Interpretation¶

Potential Use Case¶

	Age	Attrition	BusinessTravel	DailyRate	Department	DistanceFromHome	Education	EducationField	EmployeeCount	EmployeeNumber	EnvironmentSatisfaction	Gender	HourlyRate	JobInvolvement	JobLevel	JobRole	JobSatisfaction	MaritalStatus	MonthlyIncome	MonthlyRate	NumCompaniesWorked	Over18	OverTime	PercentSalaryHike	PerformanceRating	RelationshipSatisfaction	StandardHours	StockOptionLevel	TotalWorkingYears	TrainingTimesLastYear	WorkLifeBalance	YearsAtCompany	YearsInCurrentRole	YearsSinceLastPromotion	YearsWithCurrManager
0	41	Yes	Travel_Rarely	1102	Sales	1	2	Life Sciences	1	1	2	Female	94	3	2	Sales Executive	4	Single	5993	19479	8	Y	Yes	11	3	1	80	0	8	0	1	6	4	0	5
1	49	No	Travel_Frequently	279	Research & Development	8	1	Life Sciences	1	2	3	Male	61	2	2	Research Scientist	2	Married	5130	24907	1	Y	No	23	4	4	80	1	10	3	3	10	7	1	7
2	37	Yes	Travel_Rarely	1373	Research & Development	2	2	Other	1	4	4	Male	92	2	1	Laboratory Technician	3	Single	2090	2396	6	Y	Yes	15	3	2	80	0	7	3	3	0	0	0	0
3	33	No	Travel_Frequently	1392	Research & Development	3	4	Life Sciences	1	5	4	Female	56	3	1	Research Scientist	3	Married	2909	23159	1	Y	Yes	11	3	3	80	0	8	3	3	8	7	3	0
4	27	No	Travel_Rarely	591	Research & Development	2	1	Medical	1	7	1	Male	40	3	1	Laboratory Technician	2	Married	3468	16632	9	Y	No	12	3	4	80	1	6	3	3	2	2	2	2